
import geopandas as pd
import contextily as ctx # Used for contextual basemaps
import matplotlib.pyplot as plt
from geocube.api.core import make_geocube # Used for rasterizing
import os
import shapely
import imageio # Used for making animated GIFs
import numpy as np
from IPython.display import Image
from osgeo import gdal # Raster operations
import zipfile
import rasterio
import rasterio.merge
import rasterio.plot
import rasterio.mask
plt.rcParams['figure.figsize'] = (20, 20)
os.listdir("input")

/usr/local/lib/python3.8/dist-packages/geopandas/_compat.py:106: UserWarning: The Shapely GEOS version (3.8.0-CAPI-1.13.1 ) is incompatible with the GEOS version PyGEOS was compiled with (3.9.0-CAPI-1.16.2). Conversions between both will be slow.
  warnings.warn(

['lds-nz-road-centrelines-topo-150k-FGDB.zip',
 'lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip',
 'statsnzpopulation-by-meshblock-2013-census-FGDB.zip',
 'statsnz2018-census-electoral-population-meshblock-2020-FGDB.zip',
 'statsnzregional-council-2021-clipped-generalised-FGDB.zip',
 'lris-lcdb-v50-land-cover-database-version-50-mainland-new-zealand-FGDB.zip']


%%time
REGC = pd.read_file("input/statsnzregional-council-2021-clipped-generalised-FGDB.zip!regional-council-2021-clipped-generalised.gdb")
AKL = REGC[REGC.REGC2021_V1_00_NAME == "Auckland Region"].copy()
# Filter out islands
AKL["geometry"] = max(AKL.geometry.explode(), key=lambda a: a.area)
# Coordinate reference system (projection)
print(AKL.crs)
# Simplify geometry to speed up clip operations
AKL = AKL.simplify(1000).buffer(1000)
ax = AKL.to_crs(epsg=3857).boundary.plot()
ax.set_title("Auckland Region clip extent")
ctx.add_basemap(ax)

epsg:2193
CPU times: user 1.79 s, sys: 979 ms, total: 2.77 s
Wall time: 14.1 s


%%time
df = pd.read_file("zip://input/lris-lcdb-v50-land-cover-database-version-50-mainland-new-zealand-FGDB.zip!lcdb-v50-land-cover-database-version-50-mainland-new-zealand.gdb")

CPU times: user 1min 43s, sys: 2.39 s, total: 1min 46s
Wall time: 1min 46s


print(df.columns)
print(df.crs)
display(df.sample(5))

Index(['Name_2018', 'Name_2012', 'Name_2008', 'Name_2001', 'Name_1996',
       'Class_2018', 'Class_2012', 'Class_2008', 'Class_2001', 'Class_1996',
       'Wetland_18', 'Wetland_12', 'Wetland_08', 'Wetland_01', 'Wetland_96',
       'Onshore_18', 'Onshore_12', 'Onshore_08', 'Onshore_01', 'Onshore_96',
       'EditAuthor', 'EditDate', 'LCDB_UID', 'geometry'],
      dtype='object')
epsg:2193


%%time
df = pd.clip(df, AKL)

CPU times: user 45.8 s, sys: 0 ns, total: 45.8 s
Wall time: 45.8 s


df.sample(5)


df.Name_2018.value_counts()

Exotic Forest                                3981
Indigenous Forest                            3673
Manuka and/or Kanuka                         2282
Broadleaved Indigenous Hardwoods             1788
Built-up Area (settlement)                   1350
High Producing Exotic Grassland              1326
Mangrove                                     1151
Urban Parkland/Open Space                    1099
Estuarine Open Water                          441
Orchard, Vineyard or Other Perennial Crop     436
Short-rotation Cropland                       362
Lake or Pond                                  326
Herbaceous Saline Vegetation                  303
Low Producing Grassland                       291
Gorse and/or Broom                            287
Forest - Harvested                            266
Sand or Gravel                                252
Deciduous Hardwoods                           201
Surface Mine or Dump                          132
Mixed Exotic Shrubland                        120
Herbaceous Freshwater Vegetation              118
Transport Infrastructure                      107
River                                          15
Gravel or Rock                                  9
Flaxland                                        9
Fernland                                        2
Matagouri or Grey Scrub                         1
Name: Name_2018, dtype: int64


def simplify_classes(code):
    if code in [1, 2, 5]:
        return 1, "Urban"
    elif code in [68,69,71]:
        return 2, "Vegetation"
    elif code in [0,20,21,22,45,46]:
        return 3, "Water"
    else:
        return 4, "Other"

summary = []
years = [1996, 2001, 2008, 2012, 2018]
for year in years:
    print(year)
    class_year = f"Class_{year}"
    df[class_year + "_simplified_code"] = df[class_year].apply(lambda c: simplify_classes(c)[0])
    df[class_year + "_simplified_name"] = df[class_year].apply(lambda c: simplify_classes(c)[1])
    summary.append(df[class_year + "_simplified_name"].value_counts())

1996
2001
2008
2012
2018


pd.GeoDataFrame(summary).plot.area()

<AxesSubplot:>


%%capture
# %%capture suppresses output
ims = []
years = [1996, 2001, 2008, 2012, 2018]
for year in years:
    ax = df.plot(column=f'Class_{year}_simplified_name', legend=True)
    ax.set_title(year)
    ax.figure.tight_layout()
    canvas = ax.figure.canvas
    canvas.draw() # draw the canvas, cache the renderer
    image = np.frombuffer(canvas.tostring_rgb(), dtype='uint8')
    image = image.reshape(canvas.get_width_height()[::-1] + (3,))
    ims.append(image)

imageio.mimsave("land_use.gif", ims, fps=1)


with open('land_use.gif','rb') as file:
    display(Image(file.read()))


cols = [f"Class_{year}_simplified_code" for year in years]
cols

['Class_1996_simplified_code',
 'Class_2001_simplified_code',
 'Class_2008_simplified_code',
 'Class_2012_simplified_code',
 'Class_2018_simplified_code']


%%time
geocube = make_geocube(
    vector_data=df,
    output_crs="epsg:2193",
    measurements=cols,
    resolution=(-100, 100),
    fill=0, # NaNs, like offshore areas, will be 0
)
geocube

CPU times: user 21 s, sys: 9.84 ms, total: 21 s
Wall time: 21 s

<xarray.Dataset>
Dimensions:                     (x: 1001, y: 1320)
Coordinates:
  * y                           (y) float64 6.002e+06 6.002e+06 ... 5.87e+06
  * x                           (x) float64 1.704e+06 1.704e+06 ... 1.804e+06
    spatial_ref                 int64 0
Data variables:
    Class_1996_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2001_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2008_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2012_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2018_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


geocube.Class_2018_simplified_code.plot()

<matplotlib.collections.QuadMesh at 0x7f4086cf5460>


for year in years:
    print(year)
    outfile = f"output/land_use_{year}.tif"
    if not os.path.isfile(outfile):
        geocube[f"Class_{year}_simplified_code"].rio.to_raster(outfile, dtype=np.byte) # Use np.byte for smaller output filesize

1996
2001
2008
2012
2018


%%time
pop2013 = pd.read_file("input/statsnzpopulation-by-meshblock-2013-census-FGDB.zip!population-by-meshblock-2013-census.gdb")

CPU times: user 9.13 s, sys: 239 µs, total: 9.13 s
Wall time: 9.12 s


%%time
pop2013 = pd.clip(pop2013, AKL)

CPU times: user 18.7 s, sys: 0 ns, total: 18.7 s
Wall time: 18.7 s


display(pop2013.sample(5))


#pop2013.Population_Count_Usual_Resident_2013.replace(0, np.nan, inplace=True)


pop2013.Population_Count_Usual_Resident_2013.plot(kind="hist", bins=200)

<AxesSubplot:ylabel='Frequency'>


%%time
pop2013_cube = make_geocube(
    vector_data=pop2013,
    measurements=["Population_Count_Usual_Resident_2013"],
    like=geocube, # Ensures dimensions match
    fill=0 # NaNs, like offshore areas, will be 0
)
pop2013_cube

CPU times: user 2.68 s, sys: 10.1 ms, total: 2.69 s
Wall time: 2.69 s

<xarray.Dataset>
Dimensions:                               (x: 1001, y: 1320)
Coordinates:
  * y                                     (y) float64 6.002e+06 ... 5.87e+06
  * x                                     (x) float64 1.704e+06 ... 1.804e+06
    spatial_ref                           int64 0
Data variables:
    Population_Count_Usual_Resident_2013  (y, x) float64 0.0 0.0 0.0 ... 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


pop2013_cube.Population_Count_Usual_Resident_2013.plot()
outfile = "output/pop2013.tif"
if not os.path.isfile(outfile):
    # byte max value is 255, and we have larger values than that here. uint16 max value is 65535, which is fine
    pop2013_cube.Population_Count_Usual_Resident_2013.rio.to_raster(outfile, dtype=np.uint16)


%%time
pop2018 = pd.read_file("input/statsnz2018-census-electoral-population-meshblock-2020-FGDB.zip!2018-census-electoral-population-meshblock-2020.gdb")

CPU times: user 10.1 s, sys: 29.9 ms, total: 10.1 s
Wall time: 10.1 s


%%time
pop2018 = pd.clip(pop2018, AKL)

CPU times: user 21.1 s, sys: 0 ns, total: 21.1 s
Wall time: 21.1 s


display(pop2018.sample(5))


pop2018.General_Electoral_Population.replace(-999, 0, inplace=True)


pop2018.General_Electoral_Population.plot(kind="hist", bins=200)

<AxesSubplot:ylabel='Frequency'>


%%time
pop2018_cube = make_geocube(
    vector_data=pop2018,
    measurements=["General_Electoral_Population"],
    like=geocube, # Ensures dimensions match
    fill=0
)
pop2018_cube

CPU times: user 2.47 s, sys: 33 µs, total: 2.47 s
Wall time: 2.47 s

<xarray.Dataset>
Dimensions:                       (x: 1001, y: 1320)
Coordinates:
  * y                             (y) float64 6.002e+06 6.002e+06 ... 5.87e+06
  * x                             (x) float64 1.704e+06 1.704e+06 ... 1.804e+06
    spatial_ref                   int64 0
Data variables:
    General_Electoral_Population  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


pop2018_cube.General_Electoral_Population.plot()
outfile = "output/pop2018.tif"
if not os.path.isfile(outfile):
    pop2018_cube.General_Electoral_Population.rio.to_raster(outfile, dtype=np.uint16)


%%time
roads = pd.read_file("input/lds-nz-road-centrelines-topo-150k-FGDB.zip!nz-road-centrelines-topo-150k.gdb")

CPU times: user 11.7 s, sys: 20 ms, total: 11.7 s
Wall time: 11.7 s


%%time
akl_roads = pd.clip(roads, AKL)

CPU times: user 22.9 s, sys: 0 ns, total: 22.9 s
Wall time: 22.9 s


# If a road has a highway number (hway_num not None), it's a highway/motorway
mway = akl_roads[~akl_roads.hway_num.isna()].copy()
mway


mway.name.value_counts().head(50).plot(kind="barh").invert_yaxis()


mway.hway_num.value_counts().head(50).plot(kind="barh").invert_yaxis()


ax = mway.to_crs(epsg=3857).plot()
ax.set_title("Auckland Region motorways")
ctx.add_basemap(ax)


%%time
mway_cube = make_geocube(
    vector_data=mway,
    measurements=["lane_count"],
    like=geocube, # Ensures dimensions match
    fill=0, # 0 works fine here, as every mway has at least one lane
)
mway_cube

CPU times: user 198 ms, sys: 19.8 ms, total: 218 ms
Wall time: 215 ms

<xarray.Dataset>
Dimensions:      (x: 1001, y: 1320)
Coordinates:
  * y            (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x            (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref  int64 0
Data variables:
    lane_count   (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


mway_cube.lane_count.plot()
outfile = "output/mway.tif"
if not os.path.isfile(outfile):
    mway_cube.lane_count.rio.to_raster(outfile, dtype=np.byte)


src_ds = gdal.Open("output/mway.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/mway_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
                     src_ds.RasterXSize, src_ds.RasterYSize, 1,
                     gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from motorway in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
mway_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(mway_dist.shape)
plt.imshow(mway_dist)
plt.title("Distance from motorways in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")

(1320, 1001)

Text(0.5, 1.0, 'Distance (m)')


%%time
cbd = pop2018[pop2018.MB2020_V2_00 == "0433501"].copy()
cbd.geometry = cbd.geometry.buffer(1000)
cbd_cube = make_geocube(
    vector_data=cbd,
    like=geocube, # Ensures dimensions match
    fill=0
)
outfile = "output/cbd.tif"
cbd_cube.General_Electoral_Population.rio.to_raster(outfile, dtype=np.byte)
    
src_ds = gdal.Open("output/cbd.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/cbd_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
                     src_ds.RasterXSize, src_ds.RasterYSize, 1,
                     gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from motorway in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
cbd_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(cbd_dist.shape)
plt.imshow(cbd_dist)
plt.title("Distance from CBD in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")

(1320, 1001)
CPU times: user 314 ms, sys: 84 µs, total: 314 ms
Wall time: 312 ms

Text(0.5, 1.0, 'Distance (m)')


bounds = AKL.total_bounds.tolist()
bounds

[1703081.9789640256, 5870396.320936217, 1804839.668875325, 6002367.198185163]


zf = zipfile.ZipFile('input/lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip')
tiles = [file for file in zf.namelist() if file.endswith(".tif")]
tiles

['EJ.tif',
 'DM.tif',
 'EL.tif',
 'DL.tif',
 'DJ.tif',
 'FK.tif',
 'DK.tif',
 'EK.tif',
 'FL.tif']


tile_datasets = [rasterio.open(f'zip://input/lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip!{tile}') for tile in tiles]
DEM, transformation = rasterio.merge.merge(tile_datasets, bounds = bounds, res = (100,100), dtype=np.int16)


cbd_dist.shape

(1320, 1001)


# Ensure same shape
DEM = DEM[:,0:cbd_dist.shape[0],0:cbd_dist.shape[1]]
DEM.shape

(1, 1320, 1001)


print(np.nanmin(DEM), np.nanmean(DEM), np.nanmax(DEM), DEM.shape)
rasterio.plot.show(np.where(DEM>=0, DEM, np.nan), cmap='terrain')

-32767 -18180.83639996367 697 (1, 1320, 1001)

<AxesSubplot:>


meta = tile_datasets[0].meta
print(meta)
meta.update({
    "dtype": "int16",
    "height": DEM.shape[1],
    "width": DEM.shape[2],
    "transform": transformation
})
print(meta)
outfile = "output/slope.tif"
if not os.path.isfile(outfile):
    with rasterio.open(outfile, "w", **meta) as dest:
       dest.write(DEM)

{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -32767.0, 'width': 3028, 'height': 8192, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(8.0, 0.0, 1679712.0,
       0.0, -8.0, 5963776.0)}
{'driver': 'GTiff', 'dtype': 'int16', 'nodata': -32767.0, 'width': 1001, 'height': 1320, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(100.0, 0.0, 1703081.9789640256,
       0.0, -100.0, 6002367.198185163)}


!ls -Ggh output

total 21M
-rw-r--r-- 1 1.3M Apr 15 17:49 cbd.tif
-rw-r--r-- 1 2.6M Apr 15 17:49 cbd_dist.tif
-rw-r--r-- 1 1.3M Apr 15 14:25 land_use_1996.tif
-rw-r--r-- 1 1.3M Apr 15 14:25 land_use_2001.tif
-rw-r--r-- 1 1.3M Apr 15 14:25 land_use_2008.tif
-rw-r--r-- 1 1.3M Apr 15 14:25 land_use_2012.tif
-rw-r--r-- 1 1.3M Apr 15 14:25 land_use_2018.tif
-rw-r--r-- 1 1.3M Apr 15 14:37 mway.tif
-rw-r--r-- 1 2.6M Apr 15 17:49 mway_dist.tif
-rw-r--r-- 1 2.6M Apr 15 14:31 pop2013.tif
-rw-r--r-- 1 2.6M Apr 15 14:32 pop2018.tif
-rw-r--r-- 1 2.6M Apr 15 18:16 slope.tif

	Meshblock	MeshblockNumber	Population_Count_Usual_Resident_2013	Population_Count_Census_Night_2013	geometry
2302	MB 0760134	0760134	270	264	POLYGON ((1769599.837 5903578.275, 1769601.923...
16272	MB 0547200	0547200	81	78	POLYGON ((1756915.336 5917118.791, 1756912.971...
15297	MB 0462816	0462816	30	30	POLYGON ((1763651.497 5917286.924, 1763655.014...
14093	MB 0350902	0350902	57	60	POLYGON ((1753772.129 5924969.208, 1753771.960...
12761	MB 0246700	0246700	192	189	POLYGON ((1746479.340 5922342.997, 1746481.359...

	MB2020_V2_00	General_Electoral_Population	Maori_Electoral_Population	GED2020_V1_00	GED2020_V1_00_NAME	GED2020_V1_00_NAME_ASCII	MED2020_V1_00	MED2020_V1_00_NAME	MED2020_V1_00_NAME_ASCII	LAND_AREA_SQ_KM	AREA_SQ_KM	Shape_Length	geometry
12834	0207400	153	6	031	Northcote	Northcote	5	Te Tai Tokerau	Te Tai Tokerau	0.037010	0.037010	973.695953	POLYGON ((1754420.268 5925977.275, 1754423.490...
1433	0459401	198	-999	049	Tāmaki	Tamaki	3	Tāmaki Makaurau	Tamaki Makaurau	0.051139	0.051139	1235.196095	POLYGON ((1762838.877 5918046.168, 1762842.296...
16306	0538005	78	-999	024	Mt Albert	Mt Albert	3	Tāmaki Makaurau	Tamaki Makaurau	0.015605	0.015605	495.219590	POLYGON ((1755220.193 5916006.061, 1755221.611...
15329	0436902	282	-999	001	Auckland Central	Auckland Central	3	Tāmaki Makaurau	Tamaki Makaurau	0.003800	0.003800	308.039308	POLYGON ((1757498.721 5919835.306, 1757508.451...
12113	0174212	72	-999	063	Whangaparāoa	Whangaparaoa	5	Te Tai Tokerau	Te Tai Tokerau	0.030990	0.030990	792.626093	POLYGON ((1756701.346 5944643.959, 1756714.623...

	t50_fid	name_ascii	macronated	name	hway_num	rna_sufi	lane_count	way_count	status	surface	geometry
512	100120610	KAIPARA COAST HIGHWAY	N	KAIPARA COAST HIGHWAY	16	3007739	2	None	None	sealed	LINESTRING (1732000.000 5944172.070, 1732048.5...
2933	3198057	STATE HIGHWAY 1	N	STATE HIGHWAY 1	1	3027695	2	None	None	sealed	LINESTRING (1748581.508 5968975.145, 1748558.4...
2934	3198059	STATE HIGHWAY 1	N	STATE HIGHWAY 1	1	3027695	2	None	None	sealed	LINESTRING (1748171.047 5971284.152, 1748129.9...
3320	3200754	PAERATA ROAD	N	PAERATA ROAD	22	3000260	2	None	None	sealed	LINESTRING (1767236.112 5888088.508, 1767244.3...
3324	3200792	UPPER HARBOUR MOTORWAY	N	UPPER HARBOUR MOTORWAY	18	3047073	4	None	None	sealed	LINESTRING (1747954.314 5927269.837, 1747970.0...
...	...	...	...	...	...	...	...	...	...	...	...
138240	100048291	AUCKLAND-WAIWERA MOTORWAY	N	AUCKLAND-WAIWERA MOTORWAY	1	3067966	7	None	None	sealed	LINESTRING (1755881.018 5922863.734, 1755886.4...
138301	100048432	AUCKLAND-HAMILTON MOTORWAY	N	AUCKLAND-HAMILTON MOTORWAY	1	3017109	1	None	None	sealed	LINESTRING (1765115.647 5909916.697, 1765092.7...
138337	100048532	STATE HIGHWAY 1	N	STATE HIGHWAY 1	1	3027695	4	None	None	sealed	LINESTRING (1748892.089 5949596.727, 1748892.0...
138369	100048589	PORT ALBERT ROAD	N	PORT ALBERT ROAD	16	3013274	2	None	None	sealed	LINESTRING (1734173.019 5980575.187, 1734175.6...
138680	100118365	SOUTH-WESTERN MOTORWAY	N	SOUTH-WESTERN MOTORWAY	20	3018532	4	None	None	sealed	LINESTRING (1760066.252 5908184.133, 1760043.9...

Land use¶

Population density¶

Distance from major roads¶

Distance from CBD¶

Slope¶

	Name_2018	Name_2012	Name_2008	Name_2001	Name_1996	Class_2018	Class_2012	Class_2008	Class_2001	Class_1996	...	Wetland_96	Onshore_18	Onshore_12	Onshore_08	Onshore_01	Onshore_96	EditAuthor	EditDate	LCDB_UID	geometry
282583	High Producing Exotic Grassland	High Producing Exotic Grassland	High Producing Exotic Grassland	Gorse and/or Broom	Gorse and/or Broom	40	40	40	51	51	...	no	yes	yes	yes	yes	yes	Landcare Research	2014-06-30T00:00:00	lcdb2000099787	MULTIPOLYGON (((1436194.806 5262006.720, 14362...
240856	Broadleaved Indigenous Hardwoods	Forest - Harvested	Exotic Forest	Exotic Forest	Exotic Forest	54	64	71	71	71	...	no	yes	yes	yes	yes	yes	Landcare Research	2019-12-01T00:00:00	lcdb1000407353	MULTIPOLYGON (((1869683.875 5814658.055, 18696...
190935	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	54	54	54	54	54	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb2000140213	MULTIPOLYGON (((1643078.666 5304442.323, 16430...
298903	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	52	52	52	52	52	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000080516	MULTIPOLYGON (((1795333.019 5617045.693, 17953...
138646	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	52	52	52	52	52	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb2000128420	MULTIPOLYGON (((1561655.818 5484526.615, 15616...

	Name_2018	Name_2012	Name_2008	Name_2001	Name_1996	Class_2018	Class_2012	Class_2008	Class_2001	Class_1996	...	Wetland_96	Onshore_18	Onshore_12	Onshore_08	Onshore_01	Onshore_96	EditAuthor	EditDate	LCDB_UID	geometry
290956	Indigenous Forest	Indigenous Forest	Indigenous Forest	Indigenous Forest	Indigenous Forest	69	69	69	69	69	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000164146	POLYGON ((1775884.006 5889569.522, 1775878.437...
82724	Urban Parkland/Open Space	Urban Parkland/Open Space	Urban Parkland/Open Space	Urban Parkland/Open Space	Urban Parkland/Open Space	2	2	2	2	2	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000006339	POLYGON ((1756145.229 5927664.912, 1756146.259...
385423	Forest - Harvested	Exotic Forest	Exotic Forest	Exotic Forest	Exotic Forest	64	71	71	71	71	...	no	yes	yes	yes	yes	yes	Landcare Research	2019-12-01T00:00:00	lcdb1000219696	POLYGON ((1716669.127 5951801.208, 1716650.371...
507928	Exotic Forest	Exotic Forest	Exotic Forest	Exotic Forest	Exotic Forest	71	71	71	71	71	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000220692	POLYGON ((1744258.558 5989598.246, 1744248.730...
34208	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	52	52	52	52	52	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000095213	POLYGON ((1742543.122 5967517.024, 1742571.974...